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ABSTRACT 

We present a new method for constructing three-dimensional mass maps from gravitational lensing 
shear data. We solve the lensing inversion problem using truncation of singular values (within the con- 
text of generalized least squares estimation) without a priori assumptions about the statistical nature 
of the signal. This singular value framework allows a quantitative comparison between different filter- 
ing methods: we evaluate our method beside the previously explored Wiener filter approaches. Our 
method yields near-optimal angular resolution of the lensing reconstruction and allows cluster sized 
halos to be de-blended robustly. It allows for mass reconstructions which are 2-3 orders-of-magnitude 
faster than the Wiener filter approach; in particular, we estimate that an all-sky reconstruction with 
arcminute resolution could be performed on a time-scale of hours. We find however that linear, non- 
parametric reconstructions have a fundamental limitation in the resolution achieved in the redshift 
direction. 

Subject headings: gravitational lensing — dark matter — large-scale structure of universe 



1. INTRODUCTION 

In recent years, the study of gravitational lensing has 
proved a valuable tool in advancing our understanding 
of the universe. Because deflection of light in a grav- 
itational field is a well-understood aspect of Einstein's 
General Relativity, it offers a unique method of map- 
ping the mass distribution within the universe (includ- 
ing the dark matter), which is free of any astrophysical 
bias. Though much insight can be gained from the two 
dim ensional projectio n of the matter distribution (see, 
e.g. iClowe et al.ll2006f) . a nonparametric technique that 
can map the full 3D matter distribution is in principle 
a ttainabl e. 

iTavlorl (I 200P). Hu fc Keetonl (piM hereafter HK02) 
and [Bacon fc Tavlorl (|20n3T first looked at non- 
parametric 3D mapping of a gravitational potential. 
HK02 presented a linear-algebraic method for tomo- 
graphic mapping of the matter distribution - splitting the 
sources and lenses into discrete planes in redshift. They 
found that the inversion along each line-of-sight is ill- 
conditioned, and requires regularization through Wiener 
filtering. Wiener filtering reduces reconstruction noise by 
using the expected statistical properties of the signal as 
a prior: for the present probl em, this prio r is th e non- 
linear mass power spectrum. iSimon et al.l (|2009t here- 
after STH09) made important advances to this method 
by constructing an efficient framework in which the in- 
versions for every line-of-sight are computed simultane- 
ously, allowing for greater flexibility in the type of filter 
used. They introduced two types of Wiener filters: a 
"radial Wiener filter" , based on the HK02 method, and a 
"transverse Wiener filter" , based on the Limber approx- 
imation to the 3D mass power spectrum. They showed 
that the use of a generalized form of either filter leads to 
a biased result - the filtered reconstruction of the line- 



of-sight matter distribution for a localized lensing mass 
is both shifted and spread-out in redshift. 

One issue with the Wiener filter approach is the as- 
sumption of Gaussian statistics in the reconstructed sig- 
nal. In reality, the matter distribution at relevant scales 
can be highly non-Gaussian. It is possible that the red- 
shift bias found in STH09 is not inherent to nonparamet- 
ric linear mapping, but rather a result of this deficiency 
in the Wiener filtering method. 

In this work, we develop an alternate noise-suppression 
scheme for tomographic mapping that, unlike Wiener fil- 
tering, has no dependence on assumptions about the sig- 
nal. Our goal is to explore improvements in the recon- 
struction and examine, in particular, the recovery of red- 
shift information using the different methods. We begin 
in Section [2] by discussing the tomographic weak-lensing 
model developed by HK02 and STH09 and presenting 
our estimator for the density parameter, 5. In Sections [3] 
and U we implement this method for a simple case, and 
compare the results with those of the STH09 transverse 
and radial Wiener filters. 

2. METHOD 

For tomographic weak lensing, we are concerned with 
three quantities: the complex-valued shear j(6,z), the 
real-valued convergence k(8,z), and the dimensionless 
density parameter 5(6, z). The relationship between 7 
and k is given by a convolution over all angles 9, and the 
density 6 is related to k by a line-of-sight integral over 
the lensing efficiency function, W(z,z s ). The key obser- 
vation is that in the weak lensing regime, each of these 
operations is linear: if the variables are discretized, they 
become systems of linear equations, which can in princi- 
ple be solved using standard matrix methods. 
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2.1. Linear Mapping 

To achieve this, we create a common pixel binning of 
N x by N y equally sized square pixels of angular width 
A9 X 



A9 y . Within each of the N x N y 



N xy individual 



lines of sight, we bin 7 into N s source-planes, and bin 
S into Ni lens-planes, Ni < N s . Thus we have two ID 
data vectors, which are concatenations of the line-of-sight 
vectors within each pixel: 7, of length N xy N s ] and S, 
of length N xy Ni. (Note that throughout this section, 
boldface denotes a vector quantity.) As a result of this 
binning, we can write the discretized lensing equations 
in a particularly simple form: 



7 = M~,s5 + n 7 



(1) 



where 7 is the vector of binned shear observations with 
noise given by n 7 , and d is the vector of binned density 
parameter. For details on the form of the matrix M^s, 
refer to Appendix [A"l 

The linear estimator S of the signal is found by mini- 
mizing the quantity 

X 2 = (7 - M jS 5)W£ (7 - M lS 8) (2) 

where f indicates the conjugate transpose, and A/" 77 = 
(n 7 n 7 ^) is the noise covariance of the measurement 7, 
and we assume (n 7 ) = 0. T he best l inear unbiased esti- 
mator for this case is due to lAitkenl f|1934f ) : 



S A 



(3) 



The noise properties of this estimator can be made clear 

- — 1/2 

by defining the matrix M 7 s = A/" 77 M 1 s and computing 
the singular value decomposition (SVD) M 7 s = J7EF' . 
Here U^U = V*V = I and E is the square diagonal 
matrix of singular values <7j = ordered such that 
Cj > Ci+i j i > 1. Using these properties, the Aitken 
estimator can be equivalently written 



S A = FS- 1 C/W 7 -V2 7 



(4) 



It is apparent in this expression that the presence of small 
singular values oi <C o\ can lead to extremely large di- 
agonal entries in the matrix E _1 , which in turn amplify 
the errors in the estimator 8a- This can be seen for- 
mally by expressing the noise covariance in terms of the 
components of the SVD: 

VE~ 2 V*. 
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(5) 

The columns of the matrix V are eigenvectors of Afss, 
with eigenvalues a~ 2 . When many small singular values 
are present, the noise will dominate the reconstruction, 
and it is necessary to use a more sophisticated estimator 
to recover the signal. 

2.2. SVD Filtering 

One strategy that can be used to reduce this noise is to 
add a penalty function to the \ 2 that will suppress the 
large spikes in signal. This is the Wiener filter approach 
explored by HK02 and STH09. A more direct noise- 
reduction method, which does not require knowledge of 
the statistical properties of the signal, involves approxi- 
mating the SVD in Equation [4] to remove the contribu- 
tion of the high-noise modes. We choose a cutoff value 



a cu t, and determine n such that a n > cr cut > a n +v We 
then dehne the truncated matrices U n , E„, and V n , such 
that U n (V n ) contains the first n columns of U (V), and 
E n is a diagonal matrix of the largest n singular values, 
n < Umax- To the extent that <r 2 ut <C X)i=i °f > the trun- 
cated matrices satisfy 

C/„E n y n t « E/E0 = M^ s (6) 

and the signal estimator in Equation [3] can be approxi- 
mated by the SVD estimator: 



«5 svd (n) ^ KE-^A/-" 1 / 2 



7 



(7) 



This approximation is optimal in the sense that it prefer- 
entially eliminates high-noise orthogonal components in 
S (cf. equation [S]) , leading to an estimator which is much 
more robust to noise in 7. 

SVDs are often used in the context of Principal Com- 
ponent Analysis, where the square of the singular value 
is equal to the variance described by the corresponding 
principal component. The variance can be thought of, 
roughly, as a measure of the information contributed by 
the vector to the matrix in question. It will be useful for 
us to think about SVD truncation in this way. To that 
end, we define a measure of the truncated variance for a 
given value of n: 
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En 9 

E^max 9 
1=1 °i 



(8) 



such that < t> cu t < 1. If v C ut = 0, then n = n max and 
we are using the full Aitken estimator. As v cut —> 1, we 
are increasing the amount of truncation. 
In practice, taking the SVD of the transformation ma- 

trix My*y Mjs is not entirely straightforward: the ma- 
trix is of size (N xy N s ) x (N xy Ni). With a 128 x 128- 
pixel field, 20 lens-planes, and 25 source-planes, the ma- 
trix contains 1.3 x 10 11 mostly nonzero complex entries, 
amounting to 2TB in memory (double precision). Com- 
puting the SVD for a non-sparse matrix of this size is far 
from trivial. 

We have developed a technique to speed-up this pro- 
cess, which involves decomposing the matrices A/" 77 and 
Mn/6 into tensor products, so that the full SVD can be 
determined through computing SVDs of two smaller ma- 
trices: an N s x Ni matrix, and an N xy x N xy matrix. The 
second of these individual SVDs can be approximated 
using the Fourier-space properties of the 7 — > k map- 
ping. The result is that the entire SVD estimator can be 
computed very quickly. The details of this method are 
described in Appendix [Al 

3. RESULTS 

Using the above formalism, we can now explore the 
tomographic weak lensing problem using the techniques 
of Section O For the following discussion, we will use a 
field of approximately one square degree: a 64 x 64 grid 
of 1' x 1' pixels, with 25 source redshift bins (0 < z < 2.0, 
Az = 0.08) and 20 lens redshift bins (0 < z < 2.0, 
Az = 0.1). This binning approximates the expected pho- 
tometric redshift errors of future surveys. We suppress 
edge effects by increasing the noise of all pixels within 4' 
of the field border by a factor of 10 3 , effectively deweight- 
ing the signal in these pixels (cf. STH09). The noise for 
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Fig. 1. — Ordered singular values of the matrix M 1 $. The dotted lines show the values of n such that 99%, 99.9%, and 99.99% of the 
variance is preserved. The sharp drop-off near n = 60, 000 is due to the 10~ 3 deweighting of border pixels. 



where ct 7 is the 



each redshift bin is set to = o^j \J 
intrinsic ellipticity dispersion, and iVj is the number of 
galaxies in the bin. We assume tr 7 = 0.3, and 70 galaxies 
per square arcminute, with a redshift distribution given 
by 



n(z)ocz 2 exp — (z/zq) 



n.3/2 



(9) 



with zq = 0.57. We assume a flat cosmology with h 
0.7. Q M = 0.27 and = 0.73 at the present day. 

3.1. Singular Values 

The singular values of the transformation matrix for 
this configuration are depicted in Figure [TJ The step 
pattern visible in this plot is due to the fact that the 
noise across each source plane is identical, aside from the 
4' deweighted border. It is apparent from this figure that 
the large majority of the singular values are very small: 
99.9% of the variance in the transformation is contained 
in less than 1/3 of the singular values. The large number 
of very small singular values will, therefore, dominate in 
the Aitken estimator (Equation @|, leading to the very 
noisy unfiltered results seen in HK02. 

3.2. Evaluation of the SVD Estimator 

To evaluate the performance of the SVD filter, we first 
create a field-of-view containing a single halo at redshift 
z = 0.6. One well-support ed parametrization o f halo 
shapes is the NFW profile (|Navarro et al.l H997) . We 



use the analytic form of the shear and projected den- 
sity due to an NFW p rofile, given by equations 13-18 in 
iTakada fc Jalnl (|2003l ). 

We reconstruct the density map using the SVD filter 
(Figure [2]) with the above survey parameters. We show 
the results for three different values of v cu t: 0.1, 0.01, and 
0.005. In all three cases, the halo is easily detected at its 
correct location (left panels), although as v cut decreases, 
there is more noise in the surrounding field. The right 
panels show the computed density profile along the line 
of sight for the central pixel. The peak of this curve is 
reasonably close to the correct redshift, but there is a 
significant spread in redshift, as well as a bias. As the 
level of SVD filtering (measured by i> cu t) decreases, the 
magnitude of these effects decreases, but the increased 
noise leads to spurious peaks. 

Similar plots for the transverse Wiener filter recom- 
mended by STH09 are shown in the upper panels of 
Figure using their recommended value of a = 0.05. 
The response shows a significant spread in angular space, 
and the signal is seen to be suppressed by six orders-of- 
magnitude along with a similar suppression of the noise. 
These effects worsen, in general, as the filtering level a 
increases. Mathematically it is apparent why the trans- 
verse filter performs so poorly: the small singular values 
primarily come from the line-of-sight part of the map- 
ping, and the this filter has no effect along the line-of- 
sight. 

The effect of the radial Wiener filter is shown in the 
bottom panels of Figure [3j It shares the positive aspects 
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Fig. 2. — The effect of SVD truncation on a single z = 0.6 NFW halo in the center of the field, for three different levels of filtering, left 
column: reconstructed density parameter 8(9) in the z = 0.6 lens-plane. The true matter distribution is represented by a tight "dot" in 
the center of the plot, right column: line-of-sight profile at the central pixel. The grey shaded area shows the input density parameter. 
The solid line shows the E-mode signal, while the dashed line shows the B-mode signal, n gives the number of singular values used in 
the reconstruction (out of a total n max = 81920), and v cu t gives the amount of variance cut by the truncation (Equation |8j ; the level of 
filtration decreases from the top panels to the bottom panels. The bottom panels show a case of under- filtering: for small enough v cu t, the 
noise overwhelms the signal and creates spurious peaks along the line-of-sight. 
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Fig. 3. — The effect of Wiener filtering on the same input as Figurc[2] Here we have used both transverse (top panels) and radial (bottom 
panels) Wiener filtering, both down-tuned by a = 0.05 (the value recommended by STH09). The transverse Wiener filter suppresses the 
response by several orders of magnitude; a closer view of the line-of-sight peak is shown in the inset plot. The radial Wiener filter gives 
similar angular results to the SVD filter, but takes much longer to compute. 
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of the SVD filter, having very little signal suppression or 
angular spread. However, this filter uses some priors on 
the statistical form of the signal that are not as physically 
well-motivated as those for the transverse Wiener filter. 
In contrast, the SVD filter does not make any prior as- 
sumptions about the signal. In this way, the SVD recon- 
struction can be thought of as even more non-parametric 
than the Wiener filter reconstructions. 

3.3. Comparison of Estimators 

The SVD framework laid out in Section 12.21 can be 
used to quantitatively compare the behavior of different 
estimators. A general linear estimator has the form 

Sr = Rrf (fO) 

for some matrix R. This general estimator can be ex- 
pressed in terms of the components of the unbiased esti- 
mator (Equation |4]): 

r = v r z- 1 u^-v 2 . (11) 

Here the matrices S, U and A/" 77 are defined as in Equa- 
tion |4l and we have defined the matrix 

V R ee RM^ 2 UH (12) 

The rows of the matrix Tr^U^M^^ 2 provide a conve- 
nient basis in which to work: they are the weighted prin- 
cipal components of the shear, ordered with decreasing 
signal to noise. The norm of the i th column of Vr mea- 
sures the contribution of the i th mode to the reconstruc- 
tion of S. For the unfiltered estimator, Vr = V and 
all the norms are unity. This leads to a very intuitive 
comparison between different filtering schemes. Figure [4] 
compares the column- norms of Vr for the SVD filter with 
those of the radial and transverse Wiener filters. 

The steps visible in the plot originate the same way as 
the steps in Figure[TJ the flatness of each step comes from 
the assumption of uniform noise in each source plane. 
This plot shows the tradeoff between noise and bias. The 
flat line at norm=10° represents a noisy but unbiased es- 
timator. Any departure from this will impose a bias, but 
can increase signal-to-noise. There are two important 
observations from this figure. First, because each step 
on the plot is relatively flat for the SVD filter and radial 
Wiener filter, we don't expect much bias within each lens 
plane. The transverse filter, on the other hand, has fluc- 
tuations at the 10% level within each step (visible in the 
inset of Figure , which will lead to a noticeable bias 
within each lens plane, resulting in the degraded angular 
resolution of the reconstruction seen in Figure[3] Second, 
the transverse Wiener filter deweights even the highest 
signal-to-noise modes by many orders of magnitude, re- 
sulting in the signal suppression seen in Figure El The 
SVD filter and radial Wiener filter, on the other hand, 
have weights near unity for the highest signal-to-noise 
modes. These two observations show why the SVD fil- 
ter and radial Wiener filter are the more successful noise 
reduction techniques for the present problem. 

3.4. Noise Properties of Line-of-Sight Modes 

As seen in equation [SJ the columns of V provide a nat- 
ural orthogonal basis in which to express the signal <5. It 
should be emphasized that this eigenbasis is valid for any 



linear filtering scheme: the untruncated SVD is simply 
an equivalent re-expression of the original transforma- 
tion. Examining the characteristics of these eigenmodes 
can yield insight regardless of the filtering method used. 

The radial components of the first four eigenmodes are 
plotted in figure [5l Each is labeled by its normalized 
noise level, rii ee ((Ji/ai)^ 1 . The total number of modes 
will be equal to the number of output redshift bins; here, 
for clarity, we've used 80 equally-spaced bins out to red- 
shift 2.0. As the resolution is lessened, the overall shape 
and relative noise level of the lower-order modes is main- 
tained. These radial modes are analogous to angular 
Fourier modes, and are related to the signal-to-noise KL 
modes discussed in HK02. It is clear from this plot that 
any linear, non-parametric estimator will be fundamen- 
tally limited in its redshift resolution: the noise level of 
the i th mode approximately scales as 

ni oc i 2 (13) 

The signal-to-noise level for any particular halo will de- 
pend on its mass and redshift. The magnitude of the sig- 
nal scales linearly with mass (see discussion in STH09), 
but the redshift dependence is more complicated: it is 
affected by the lensing efficiency function, which de- 
pends on the redshift of the lensed galaxies. Using the 
above survey parameters, with an NFW halo of mass 
= 10 15 M© and redshift z = 0.6, the signal-to-noise 
ratio of the central pixel for the fundamental radial mode 
is ~ 5.9, consistent with the results for Wiener filtered 
reconstructions of singular isothermal halos explored in 
STH09. This means that for even the largest halos, with 
a very deep survey, only the first few modes will con- 
tribute significantly to the reconstructed halo. Adding 
higher-order modes can in theory provide redshift infor- 
mation, but at the cost of increasingly high noise con- 
tamination. This is a general result which will apply to 
all nonparametric linear reconstruction algorithms. 

This lack of information in the redshift direction leads 
directly to an inability to accurately determine halo 
masses: the lensing equations relate observed shear 7 
to density parameter <5, which is related to mass in a 
redshift-dependent way. This is a fundamental limita- 
tion on the ability of linear nonparametric methods to 
determine halo masses from shear data. Indeed, even 
moving to fully parametric models, line-of-sig ht effects 
can lead to halo mass errors o f 20% or more (|Hoekstral 
l200l:lde Putter fc Whitdl2005h . 

3.5. Reconstruction of a Realistic Field 

To compare the performance of the three filtering 
methods for a realistic field, we create a 4 square de- 
gree field with approximately 20 halos between masses 
of 2 x 10 14 and 8 x 10 14 M Q with a mass distribution 
appro ximating the cluster mass function of iRines et al.l 
(2007), and a redshift distribution given by Equation^ 
adding a hard cutoff at z = 1.0. These parameters are 
chosen to approximate the true distribution of observable 
halos in a field this size. The results of the reconstruction 
are shown in Figure |6] 

The red circles are the locations of the input halos, not 
the result of some halo-detection algorithm. However, it 
is clear that, for at least most of the mass range, we are 
able to produce a map for which any reasonable detection 
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Fig. 4. — Contribution of each shear mode to the reconstruction for three different filters. The dotted line at 10° represents the unfiltered 
result. Each filtering method leads to a different weighting of the shear modes. The SVD filter, by design, completely removes higher-order 
modes beyond a given cutoff, while the Wiener filter deweights modes in a more gradual fashion. Note that the transverse Weiner filter 
deweights all modes by up to seven orders of magnitude; it has been scaled by a factor of 10 4 for this plot. The inset plot shows a closeup 
of the fluctu ations within each "step" of the transverse filter. These fluctuations lead to angular spread in the response (see discussion in 
Section 13.31 1 



algorithm should detect the halos in the correct locations. 
A few of the lower mass halos would certainly be missed 
though, since they are not significantly different from the 
noise peaks in the image. 

In practice, one may vary the parameter u cut as in Fig- 
ure [5] to trade-off robustness of detecting peaks with res- 
olution in angle and in redshift. As shown in Section f3.3[ 
we expect filtering to introduce very little bias in angular 
resolution, so large values of v cut lead to the most robust 
angular results. On the other hand, as shown in Sec- 
tion 13.41 filtering introduces an extreme bias along the 
line-of-sight. The effects of this bias can be seen qualita- 
tively in the right column of Figure [2] Optimal redshift 
resolution requires choosing a filtering level which bal- 
ances the effects of noise and bias, and may require some 
form of bias correction. In future work, we will explore 
in detail the ways in which the SVD method allows for a 
near optimal reconstruction of projected mass maps and 
halo redshifts from data on galaxy shapes and photomet- 
ric redshifts. 

3.6. Scalability 

As we look forward to future surveys, it becomes im- 
portant to consider methods that will scale upward with 
increasing survey volumes. Present weak lensing surveys 
cover field s on the order of a few square degrees (e.g. 
COSMOS. iMassev et all 120071) . Future surveys will in- 



crease the field size e xponentially: up to ~20, 000 square 
degre es for LSST (jLSST Science Collaborations et al.l 
2009). Though the fiat-sky approximation used in this 
work is not appropriate for such large survey areas, the 
weak lensing formalism can be modified t o account for 
spherical geometry (see, e.g. lHeavensll2003l) . 

The main computational cost for both SVD and 
Wiener filtering is the Fast Fourier Transform (FFT) 
required to implement the mapping from 7 to n. For 
an N x N pixel field, the FFT algorithm performs in 
0[N log N] in each dimension, meaning that the 2D FFT 
takes 0{{N log A) 2 ] w 0[N 2 ]. The Wiener filter method, 
however, requires the inversion of a very large matrix us- 
ing, for example, a conjugate-gradient method. The ex- 
act number of iterations depends highly on the condition 
number of the matrix to be inverted; STH09 finds that 
up to 150 iterations are required for this problem. We 
find that each iteration takes over 3 times longer than the 
entire SVD reconstruction. The net result is that both 
algorithms scale nearly linearly with the area of the field 
(for constant pixel scale), though the SVD estimator is 
computed up to 500 times faster than the Wiener filter. 

Extrapolating this scaling, the appropriately scaled 
SVD filter will allow reconstruction of the entire ^20, 000 
square-degree LSST field in a few hours on a single work- 
station, given enough memory. On the same computer, 
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Fig. 5. — left panel: The radial components of the first four columns of the matrix V (see section \2.l$ . This is calculated for 100 equally 
spaced redshift bins (0 < z < 2.5) in *y, and 80 bins (0 < z < 2.0) in S These orthogonal eigenmodes are analogous to radial Fourier modes. 
Each is labeled by its relative noise level, n, = (cj/cri) - . right panel: The singular values Cj associated with the 80 radial eigenmodes. 



TABLE 1 

Masses and redshifts of halos in Figure [6] 



6 X e y z M/M e 



A 


37.5 


44.9 


0.60 


7.2 x 10 14 


B 


67.1 


70.5 


0.47 


6 x 10 14 


C 


46.9 


106.5 


0.63 


5.5 x 10 14 


D 


108.9 


94.3 


0.63 


5.4 x 10 14 


E 


97.9 


63.6 


0.39 


4.9 x 10 14 


F 


102.4 


84.8 


0.70 


3.8 x 10 14 


G 


77.0 


49.6 


0.58 


3.2 x 10 14 


H 


52.0 


48.5 


0.36 


3.2 x 10 14 


1 


72.6 


45.6 


0.78 


2.9 x 10 14 
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the Wiencr-filtcr method would take over a month, de- 
pending on the amount and type of filtering and assum- 
ing that the required number of iterations stays constant 
with increasing field size. For the SVD-filtered recon- 
struction of this large field, the real challenge will not be 
computational time, but memory constraints: the com- 
plex shear vector itself for such a field will require ~30 
GB of memory, with the entire algorithm consuming ap- 
proximately three times this. The memory requirements 
for the Wiener filter will be comparable. This is within 
reach of current high-end workstations as well as shared- 
memory parallel clusters. 

4. CONCLUSION 

We have presented a new method for producing to- 
mographic maps of dark matter through weak lensing, 
using truncation of singular values. We have tested 
and compared our method to the Wiener filter based 
method of STH09, which is the first three-dimensional 
mass mapping approach that is applicable to large area 
surveys. Our reconstruction shares many of the aspects 
of the Wiener filter reconstruction, in the sense that 
it massively reduces the noise inherent in the problem. 
Our SVD method may be considered even more non- 
parametric than the Wiener filter method, since it does 
not rely on any a priori assumptions of the statistical 



properties of the signal: all of the noise reduction is de- 
rived from the observed noise properties of the data. 

The SVD framework allows a unique quantitative com- 
parison between the different filtering methods and filter- 
ing strengths. Using the coefficients of the weighted prin- 
cipal components contained in the SVD, we have com- 
pared the three filtering methods, and have found that 
the radial Wiener filter of HK02 and SVD filter of this 
work are less-biased noise reduction techniques than the 
transverse Wiener filter of STH09. These authors have 
recently implemented the radial Weiner filter and obtain 
results consistent with our findings (P. Simon and A. 
Taylor, private communication). 

The angular resolution of the SVD-reconstructed mass 
maps seems to be significantly better than that of the 
transverse Wiener filter method, the method chosen in 
the STH09 analysis. This allows for more robust sep- 
aration of pairs of halos into two separate halos rather 
than blurring them into a single mass peak. We discuss 
how our reconstruction method provides a scheme for op- 
timizing the 3D reconstruction of projected mass maps 
by balancing the goals of robustness of detecting specific 
structures and improved redshift resolution. 

The SVD method can compute the three-dimensional 
mass maps rapidly provided sufficient computational 
memory is available. This allows for the possibility of 
solving the full-sky tomographic lensing inversion on the 
scale of hours, rather than months, which makes it read- 
ily applicable to upcoming surveys. 

On the other hand, the redshift resolution with the 
SVD method is not significantly better than that of ei- 
ther Wiener filter method. This was a problem identi- 
fied by STH09, and unfortunately the SVD method does 
not significantly improve the situation. Our analysis of 
the noise characteristics of radial modes indicates that 
linear, non-parametric reconstruction methods are fun- 
damentally limited in this regard. 
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ful suggestions. Support for this research was provided 
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Fig. 6. — Reconstruction of an artificial shear field with the SVD filter (top panels), Transverse Wiener filter (middle panels), and Radial 
Wiener filter (bottom panels). The left column shows the projected density reconstruction across the field using each method, all smoothed 
with a 1-pixel wide Gaussian filter. Red circles indicate the true locations of the input halos. The right column shows the line-of-sight 
distributions of the twelve most massive NFW halos, labeled A-L. The masses and redshifts of the halos are listed in Table [TJ The signal 
suppression of the transverse Wiener filter seen in Figure|3]is apparent in the color-bar scaling of the middle pan els. The anomalous results 
seen in halo K are due to its proximity to the deweighted border. As suggested by the discussion in Section 13.41 none of the three methods 
succeed in recovering precise redshifts of the halos. 



APPENDIX 

EFFICIENT IMPLEMENTATION OF THE SVD ESTIMATOR 

* — — - 1/2 

As noted in Section taking the SVD of the transformation matrix M 1 s = A/" 77 M 1 $ is not trivial for large fields. 
This appendix will first give a rough outline of the form of M^s, then describe our tensor decomposition method which 
ena bles quick calculation of the sing ular value decomposition. For a more thorough review of the lensing results, see 
e.g. iBartelmann fc Schneiderl (|2001l ). 

Our goal is to speed the computation of the SVD by writing M 7 g as a tensor product A® B. Here "(g)" is the 
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Kronecker product, defined such that, if A is a matrix of size n x m, B is a matrix of arbitrary size, 



A®B 



( A n B A 12 B ••• A lm B\ 
A 2l B A 22 B ••• A 2m B 



(Al) 



\ A n iB A n2 B ■ ■ ■ A nm B J 

In this case, the singular value decomposition A <£> B = U ab^abV'ab satisfies 

Uab = U a ®U b 

Vab = V a ®V b (A2) 

where Ua^aV\ is the SVD of A, and Ub^bV^ is the SVD of B. Decomposing M 7 s in this way can greatly speed the 
SVD computation. 

Angular and Line-of-Sight Transformations 

The transformation from shear to density, encoded in M^s, consists of two steps: an angular integral relating shear 
7 to convergence k, and a line-of-sight integral relating the convergence k to the density contrast S. 
The relationship between 7 and k is a convolution over all angular scales, 

7(0, z„) = 7l + ^7 2 = J dV V(0' - 8)k{8\ z s ), (A3) 

where T>(9) is the Kaiser- Squires kernel ()Kaiser fc Squ ires 1993). This has a particularly simple form in Fourier space: 

K£, Zs ) = *±±^Lfi(l, Zs ). (A4) 

where 7 and k are the Fourier transforms of 7 and k and £ = (£1, £ 2 ) is the angular wavenumber. 
The relationship between n and S is an integral along each line of sight: 

k(0,z s )= f dz W(z,z s )6(6,z) (A5) 
Jo 

where W(z, z s ) is the lensing efficiency function at redshift z for a source located at redshift z s (refer to STH09 for 
the form of this function). 

Upon discretization of the quantities 7, k, and 6 (described in Section COj) . the integrals in Eq uations I A3llA"5l become 
matrix operations. The relationship between the data vectors 7 and k can be written 

7 = [P~fn <£> l s ]« + n~f (A6) 

where l s is the N s x iV s identity matrix and P 7K is the matrix representing the linear transformation in Equations lA3t - 
IA4I The quantity [P 7K l s ] simply denotes that P 7K operates on each of the N s source-planes represented within the 
vector k. Similarly, the relationship between the vectors k and S can be written 

k = [l xy <g) Q kS ]8 (A7) 

where l xy is the N xy x 7V xy identity matrix, and the tensor product signifies that the operator Q K $ operates on each 
of the N xy lines-of-sight in 6. Q K $ is the N s x Ni matrix which represents the discretized version of equation IA51 
Combining these representations allows us to decompose the matrix My$ in Equation [T] into a tensor product: 

M 1& = P 7K ® Q kS . (A8) 

Tensor Decomposition of the Transformation 

We now make an approximation that the noise covariance A/" 77 can be written as a tensor product between its angular 
part Afp and its line of sight part ATq: 

N 11 =N'p®N'q. (A9) 

Because shear measurement error comes primarily from shot noise, this approximation is equivalent to the statement 
that source galaxies are drawn from a single redshift distribution, with a different normalization along each line-of-sight. 
For realistic data, this approximation will break down as the size of the pixels becomes very small. We will assume 
here for simplicity that the noise covariance is diagonal, but the following results can be generalized for non-diagonal 
noise. Using this noise covariance approximation, we can compute the SVDs of the components of M^g: 

U P ^pVl=Up l,2 P 1K 

-1/2 



UqVqV^ =AT q l/ 'Q kS (A10) 
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In practice the SVD of the matrix P 7K need not be computed explicitly. P 7K encodes the discrete linear operation 
expressed by Equations I A3II A4t as pointed out by STH09, in the large-field limit P 7K can be equivalently computed in 
either real or Fourier space. Thus to operate with P 7K on a shear vector, we first take the 2D Fast Fourier Transform 
(FFT) of each source-plane, multiply by the kernel (li + i^a) / (^1 — H2), then take the inverse FFT of the result. This is 
orders-of-magnitude faster than a discrete implementation of the real-space convolution. Furthermore, the conjugate 
transpose of this operation can be computed by transforming i — > —£*, so that 

P^Py, = I (AH) 
and we see that P 7K is unitary in the wide-field limit. This fact, along with the tensor product properties of the SVD, 
allows us to write M^s = UY0V* where 

Unl xy <g) U Q 

V^P JK ®V% (A12) 

The only explicit SVD we need to calculate is that of Nq Q k s, which is trivial in cases of interest. The two 

approximations we have made are the applicability of the Fourier-space form of the 7 — > k mapping (Eqn. IA4|) . and 
the tensor decomposition of the noise covariance (Eqn. IA9|) . 
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